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Abstract 

Transitive consistency is an intrinsic property for collections of linear invertible transformations between Euclidean coordinate 
frames. In practice, when the transformations are estimated from data, this property is lacking. This work addresses the problem 
of synchronizing transformations that are not transitively consistent. Once the transformations have been synchronized, they 
satisfy the transitive consistency condition - a transformation from frame A to frame C is equal to the composite transformation 
of first transforming A to B and then transforming B to C. The coordinate frames correspond to nodes in a graph and the 
transformations correspond to edges in the same graph. Two direct or centralized synchronization methods are presented for 
different graph topologies; the first one for quasi-strongly connected graphs, and the second one for connected graphs. As an 
extension of the second method, an iterative Gauss-Newton method is presented, which is later adapted to the case of affine 
and Euclidean transformations. Two distributed synchronization methods are also presented for orthogonal matrices, which 
can be seen as distributed versions of the two direct or centralized methods; they are similar in nature to standard consensus 
protocols used for distributed averaging. When the transformations are orthogonal matrices, a bound on the optimality gap 
can be computed. Simulations show that the gap is almost tight, even for noise large in magnitude. This work also contributes 
on a theoretical level by providing linear algebraic relationships for transitively consistent transformations. One of the benefits 
of the proposed methods is their simplicity - basic linear algebraic methods are used, e.g., the Singular Value Decomposition 
(SVD). For a wide range of parameter settings, the methods are numerically validated. 
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1 Introduction 

Collections of linear invertible transformations between 
Euclidean coordinate systems must be transitively con¬ 
sistent. In practice however, when the transformations 
are estimated from data, this condition does not hold. 
This issue is present in the 3D localization problem, 
where transformations are rigid and estimated from e.g., 
camera measurements; in the multiple images registra¬ 
tion problem where the transformations are affine (or 
linear by using homogeneous coordinates); in the gen¬ 
eralized Procrustes problem where scales, rotations and 
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translations are calculated from multiple point clouds. In 
order to resolve the issue, the estimated transformations 
need to be synchronized in the sense of finding transi¬ 
tively consistent transformations close to the estimated 
ones. 


1.1 Problem 

This work addresses the problem of synchronizing lin¬ 
ear invertible transformations or matrices between Eu¬ 
clidean coordinate systems or frames. More precisely, 
given a collection {Gij} of matrices in GL{d, K.), another 
collection {GU} of matrices in GL{d,R) is constructed 
such that 

G*jG*f, = G*fc, for all i,j, k, 


( 1 ) 





where G*j is “close” to Gij for all i,j- By satisfying (1), 
the collection {G*j} is said to be transitively consistent. 

1.2 Background 

There are many applications for the proposed methods. 
One such application is the 3D localization problem in 
camera networks [1] where a network of cameras are ob¬ 
serving a scene and epipolar geometry is used to cal¬ 
culate/measure Gij transformations between (i,j)-pairs 
of cameras. If the cameras are fully calibrated these 
transformations are Euclidean, otherwise they conld be 
e.g., affine (or linear by using homogeneous coordinates). 
Since the transformations are calculated from measure¬ 
ments, they do not satisfy (1) in general. Hence our pro¬ 
posed methods can be used to synchronize the matrices. 
For the 3D localization problem, we do not have to limit 
ourselves to the case of cameras and epipolar geome¬ 
try. The transformations could be calculated in a setting 
where the geometry of the scene is known. In the case of 
known point features, the perspective-n-point problem 
can be solved in order to get estimates of the relative 
transformations [2]. 

Another important problem is image registration, which 
has attracted much attention in the medical imaging 
community. The number of applications is vast, rang¬ 
ing from surgery planning to longitudinal studies. To 
register a (moving) image with another (fixed) image is 
to transform the former into the the latter in such way 
that they fit in the “best” way. For that, optimization 
methods are used to calculate a transformation which 
minimizes a suitable objective. Registration of multiple 
images poses a greater challenge. There are several ap¬ 
proaches in the literature. For example: Finding a path 
of pairwise transformations, which contains all images 
[3]; aligning images with a reference frame [4]; image 
congealing, where variability along known axes of vari¬ 
ation is removed in an iterative manner [4]; consider¬ 
ing a minimum description length (MDL) approach of a 
statistical shape model built from the correspondences 
given due to groupwise image registration [5]; Bayesian 
formulations and Expected Maximization (EM) [6]. 

Another way to solve the (affine) multiple images regis¬ 
tration problem is to use the transitive consistency cri¬ 
terion (I) [7]. Let the Gij correspond to the affine trans¬ 
formations calculated from pairwise registrations, then 
our method can be used to create transitively consistent 
Gij transformations. Registration methods using transi¬ 
tive consistency have also been proposed for deformable 
transformations [8,9]. 

A related problem to the one posed in this paper is the 
problem of calculating the “best” translations, rotations 
and scales between pairs of point clouds. If only one pair 
is considered the problem is referred to as the Procrustes 
Problem [10]. This problem can be solved by means of 


singular value decomposition or eigenvalue decomposi¬ 
tion [11,12,13], or in the case case of 3D transformations, 
by a quaternion-based approach [14,15]. The problem 
restricted to 3D is referred to as the absolute orientation 
problem [13,14]. In the general setting, when n point 
clouds are considered, the problem is referred to as the 
Generalized Procrustes Problem [10]. In order to solve 
this problem, iterative methods are often used; when the 
dimension is two or three, direct methods have recently 
been proposed [16]. Our previous work in [17] has tackled 
the Generalized Procrustes Problem using an approach 
based on transitive consistency. The present paper will 
extend and generalize these ideas as well as describe 
many theoretical properties of the generalizations. 

Our methods can be used for solving the Generalized 
Procrustes Problem in the following way: Between each 
pair of point clouds a Gij transformation is calculated 
using any standard technique [11,12,13], then our meth¬ 
ods are used to improve the pairwise transformations by 
calculating transitively consistent transformations. 

In the special case when the Gij are orthogonal matrices, 
Singer et al. have presented methods for the optimiza¬ 
tion of transitive consistency [18,19,20,21]. These works 
were later adapted by Pachauri et al. to the special case 
when the Gij are permutation matrices [22]. In the latter 
work, a relaxation of the original problem is considered - 
in the original problem the transformations shall be or¬ 
thogonal matrices - and then permutation matrices are 
obtained by means of projection from the solution of the 
relaxed problem. The method presented by Singer et al. 
is said to be a synchronization method for minimization 
of transitive consistency errors - a formalism adopted in 
this work. 

1.3 Methods and results 

The approach in this work share similarities with the ap¬ 
proaches of Singer et al. and Pachanri et al; it continues 
along the lines of the the recently proposed methods in 
[17,23]. 

In [17,23] a so called Z-matrix is constructed from the 
Gij matrices. If the index set for the (available) transfor¬ 
mations has a certain property, transitively consistent 
transformations can be obtained by a method where the 
Singular Value Decomposition (SVD) is calculated for 
Z. The property that must be fulfilled for the index set 
{(b j)} = is that it is the edge set of a quasi-strongly 
connected (QSC) directed graph (see Definition 2). In 
the Z-matrix approach, a set of linear algebraic equa¬ 
tions are formulated - equations which shall be satisfied 
for the case of transitively consistent transformations. 
When the transformations are not transitively consis¬ 
tent, the problem is solved in the sense of least squares 
minimization. 
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As we will show in this work, the Z-matrix appears in 
the construction of a Hessian matrix H for a quadratic 
convex function of the Gij, and under certain conditions 
it holds that H = Z + From the SVD of the Hessian 
matrix H, transitively consistent transformations can be 
calculated in the same manner as for the Z-matrix. The 
justihcation for using the i7-matrix stems from the fact 
that it is the Hessian matrix of the objective function 
in a relevant optimization problem. The justification of 
using the Z-matrix stems purely from the linear alge¬ 
braic constraints that should be satisfied for transitively 
consistent transformations. 

The Z-matrix method and the i/-matrix method are 
both direct methods, i.e., the solution is found at once. 
As an extension we also propose an iterative Gauss- 
Newton method, which uses the solution from the H- 
matrix method as initialization. For orthogonal matrices 
one can prove that this iterative scheme cannot decrease 
the objective function at all. The Gauss-Newton method 
is also adapted to the cases of affine and Euclidean trans¬ 
formations. In this case - as opposed to the result for 
orthogonal matrices - significant improvement over the 
Z-matrix method and the iJ-matrix method can be seen 
in numerical simulations. 

Many properties of the Z-matrix and the i/-matrix are 
proved in this work. For example it is shown that transi¬ 
tive consistency in the case of connected graphs is equiv¬ 
alent to the condition that the nullspace of H has dimen¬ 
sion d. Furthermore, the transitively consistent transfor¬ 
mations can be obtained as the dx d blocks in a matrix, 
the columns of which span the nullspace of H. For the 
Z-matrix only a weaker condition is formulated; if the 
graph is QSC and the transformations are transitively 
consistent, the transformations can be obtained as the 
d X d blocks in a matrix, the columns of which span the 
nullspace of Z. 

Now, in most aspects the iJ-matrix approach seems to be 
superior to the Z-matrix approach. However, one large 
benefit of using the Z-matrix over the i7-matrix is that 
it can be used in a distributed algorithm when the com¬ 
munication graph is directed. 

In a later part of the paper, two distributed methods 
are introduced for the case of orthogonal Gij transfor¬ 
mations. The hrst method is using the Z-matrix un¬ 
der the assumption that the communication graph is 
directed and QSC. The other method is using the H- 
matrix under the assumption that the communication 
graph is symmetric. The performance of the two meth¬ 
ods are almost the same in numerical experiments. The 
distributed methods are similar in structure to linear 
consensus protocols [24,25,26,27,28]. Key differences to 
those approaches is that the states here are matrices in¬ 
stead of vectors, and the states combined converge to a 
d-dimensional linear subspace instead of the consensus 
set. 


The distributed iterative methods are introduced mainly 
with communication between agents in mind, e.g., in 
networks of robots with limited communication range, 
where the robots only communicate with their neigh¬ 
bors (directly or indirectly). However, a further scenario 
of the distributed methods is parallelisation in order to 
better deal with the computational burden in the case 
of very large problem instances. 


When it is known that the transitively consistent trans¬ 
formations are orthogonal matrices, i.e., elements of 
0{d) = {R-. R & M.'^^'^,R^R = /}, a method is pro¬ 
vided for calculating an upper bound on the optimality 
gap. In the case when the Gij are also orthogonal, sim¬ 
ulations show that this gap is almost tight. As an exam¬ 
ple, for n = 100 coordinate systems, dimension d = 3, 
and randomly generated Gij matrices in 0(3), the gap 
is smaller than a tenth of a percent in average. There 
are (and will be even more in the future) applications 
where large networks of cameras, robots, satellites or 
unmanned vehicles, need to synchronize their pairwise 
relative rotations. In such applications methods that are 
near optimal and run almost in real time are of utmost 
importance to have. 


1-4 Outline 


The paper proceeds as follows. In Section 2, graphs and 
properties thereof are introduced, followed by the intro¬ 
duction of the Gij transformations and their connections 
to the graphs. We have chosen to incorporate graphs in 
the very definition of transitive consistency. Section 3 ad¬ 
dresses linear invertible transformations. In Section 3.1, 
the Z-matrix is introduced, followed by a collection of 
results and a least squares method. In Section 3.4, the 
id-matrix is introduced; in the same manner as in Sec¬ 
tion 3.1, a collection of results is provided in conjunc¬ 
tion with an algorithm. In Section 3.7 a Gauss-Newton 
method is presented, where the matrices obtained from 
the id-matrix method are used as initialization. Section 4 
consider the special case of orthogonal matrices, i.e., el¬ 
ements of 0{d). The section starts with some bounds 
on the optimality gap, and continues in Section 4.1 with 
the introduction of distributed algorithms. The reader 
interested in the distributed methods can go directly to 
this section and consult the earlier sections only for ref¬ 
erence. Section 4.2 is a small detour, where a gradient 
flow method is presented for orthogonal matrices. This 
method is employed as a baseline method, used for com¬ 
parison in some of the simulations in Section 5 - the 
section where the proposed methods are thoroughly nu¬ 
merically evaluated. 
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2 Preliminaries 

2.1 Directed Graphs 

Let G = {V,£) be a directed graph, where V = 
{1,2,..., n} is the node set and f C V x V is the edge 
set. The set M is dehned by 

M = {j : (i,j) G £}. 


The adjacency matrix A = [Aij] for the graph G is de¬ 
fined by 

_ri 

else. 

The graph Laplacian matrix is defined by 
L = diag(Al„) - A, 

where 1„ G M" is a vector with all entries equal to 1. 
In order to emphasize that the adjacency matrix A, the 
Laplacian matrix L and the A/i sets depend on the graph 
G, we may write A(G), L(G) and Ni{G) respectively. For 
simplicity however, we mostly omit this notation and 
simply write A, L, and A/i. 

Definition 1 (connected graph, undirected path) 

The directed graph G is connected if there is an undirected 
path from any node in the graph to any other node. An 
undirected path is defined as a (finite) sequence of unique 
nodes such that for any pair (i,j) of consecutive nodes in 
the sequence it holds that 

((z,j) G 8) oriij,i) G 8). 

Definition 2 (quasi-strongly connected graph, center, 
directed path) 

The directed graph G is quasi-strongly connected (QSC) 
if it contains a center. A center is a node in the graph to 
which there is a directed path from any other node in the 
graph. A directed path is defined as a (finite) sequence of 
unique nodes such that any pair of consecutive nodes in 
the sequence comprises an edge in 8. 

Definition 3 (symmetric graph) 

The directed graph G = (V, 8) is symmetric if 

{{i,j) € 8) ^ ((j,z) G 8) for all {i,j) G V x V. 

Given a graph G = iyi8), the graph G = (y,8) is the 
graph constructed by reversing the direction of the edges 
in 8, i.e., (z, j) G f if and only if (j, z) G 8. It is easy to 
see that 

A{G) = {A{G)f and L{G) = dia.g{{A{G)f U) - A{Gf. 


2.2 Transformations 

Given a directed graph G = (V, 8), let there be a collec¬ 
tion of matrices where Gij G GL{d,'R) for 

all {i,j)£ 8. Let n = |V|. The Gij are not necessarily 
transitively consistent in that 

Gik GijGjk 

may hold if (z, j), (j, k) and (z, k) are elements of 8. 

In the methods to be defined, the goal is to find a tran¬ 
sitively consistent collection {G*j}(^i j'f^£ of matrices in 
GL{d,R), such that for all (z, j) G 8, G*j is close to Gij 
in some appropriate sense. Notation-wise, G*j is simply 
(a name of) a matrix. This notation should not be mixed 
up with the conjugate transpose - in this paper, all ma¬ 
trices considered are real and the conjugate transpose 
will not be used. 

Definition 4 (transitive consistency) 

(1) The matrices in the collection {G*j}(j^j)gVxv of ma¬ 
trices in GL{d, R) are transitively consistent for the 
complete graph if 

^ik — 

for all i,j and k. 

(2) Given a graph G = {V,8), the matrices in the 

collection {G*j}(ij)^£ of matrices in GL{d,R.) 
are transitively consistent for G if there is a col¬ 
lection {GL}(ij)gvxv {G*j}{i,j)ee such that 

j)gVxv is transitively consistent for the 
complete graph. 

If it is apparent by the context, sometimes we will be less 
strict and omit to mention which graph a collection of 
transformations is transitively consistent for. A sufficient 
condition for transitive consistency of the G*j matrices 
for any graph is that there is a collection {G*}igv of 
matrices in GL{d,R.) such that 

ij ^ j 

for all z, j. Lemma 6 below and the proof thereof provides 
additional important information. The result is similar 
to that in [1]. For the statement of the lemma, the fol¬ 
lowing definition is needed. 

Definition 5 Two collections {G*}igv and {G**}igv of 
matrices in GL{d, R) are equal up to transformation from 
the left, if there is Q € GL{d,R.) such that 

QG* = G*i* for all i. 
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Lemma 6 For any graph Q = {V,S) and collection 
{^ij}{i,j)€£ of matrices in GL(fi, M) that are transitively 
consistent for Q, 

(1) there is a collection {G* }iev of matrices in GL{d, K) 
such that 


G*^=G*-^G*forall{z,j)G£, (2) 

(2) all collections {G*}i^v satisfying (2) are equal up 
to transformation from the left if and only if Q is 
connected, 

(3) there is a unique collection {G*j}(i j)gVx V 3 
{^ij}transitively consistent matrices for 
the complete graph, if and only if all collections 
{G*}igv satisfying (2) are equal up to transforma¬ 
tion from the left. 

Proof: All matrices appearing in this proof, if the con¬ 
trary is not explicitly stated, are assumed to be elements 
of GL{d,R). 

(1) Since the matrices in {G*j}(j are transitively 
consistent for G, there is {Gy}(ij)gvxv 
in which the matrices are transitively consistent for the 
complete graph. Let the G* in a collection {G*}igv be 
defined by 

i — 

We shall prove that 

G*-1g* = G*-iG*y = Gb for all t,j. 

Using the fact that G^^ invertible and the fact that 
G*i = Gfl, one can show that Gff = I. Now, Gl^G*i = 
Gii = /; thus Gl~^ = G*i. But then 

G * — 1^* _ 

li — ^ij' 


(2) We know that transitive consistency of {Gb}(ij)g£ 
for Q is equivalent to the statement that there is a col¬ 
lection {G*}igv of matrices in GL{d,R) such that 

Gb = G^^G* for all (t, j) € £. 

Let Gb = Gf^Gi for alH, j G V. For any other collection 
{G**}igv of matrices in GL(d, R) such that 


Gb = G"” Gf for all (t, j) G S, 


it holds that 


G **T 

1 ^2 • • • 


— diS-g (^ 11^25 ■ • ■ ^ Qn^ 


iT 


Gt^ G*2^ ... G 


*T 


iT 


where the Qi matrices are elements of GL{d,M.). 

If: Now, if the graph is connected and at least two of the 
Qi are not equal, there is (j, k) G £ such that Qj ^ Qk. 
We know 

but since Qj ^ Qk we can calculate this entity to 

G*G*,G*-' = Q-^Qk Id, 

which is a contradiction. Q G is the identity ma¬ 

trix. 

Only if: On the other hand, if the graph is not connected 
there are two disjoint sets Vi and V 2 such that Vi U V 2 = 
V, for which there is no pair {i,j) G £ such that (i G Vi 
and j G V 2 ) or (j G Vi and i G V 2 ). Thus, the nodes 
in Vi and the corresponding edges, respective the nodes 
in V 2 and the corresponding edges, can be seen as two 
different disconnected (sub)graphs, each of them being 
connected; the G* matrices in the first graph can be 
multiplied with a matrix Qi from the left and the G* 
matrices in the second graph can be multiplied with a 
matrix Q 2 from the left, where Qi ^ Q 2 , generating a 
collection of matrices {G**}igv not equal to {G*}igv up 
to transformation from the left. 


(3) If: Any other collection {G*|y} of matrices in 
GL(fi, R) such that 

Gb = G**-^G*/ for all {i,j) G £, 

is equal to {G*} up to transformation from the left. Now, 
for any (f, j) G (V x V) — £1 it holds that 

G **—I/O** /O*—1/O —1/O/O* /O*—I/O* /O* 

i ^ j ^ ^ j 

for some matrix Q G GL{d, R). 


Only if: The approach here is similar to that in 2) above. 
Suppose for {G*gy} satisfying (2), there is another col¬ 
lection {G*|y} of matrices in GL((i, R) also satisfying 
(2), but the matrices in the two collections are not equal 
up to transformation from the left. Then it holds that 


iT 


^**T ^* 


.. G 


**T 


= diag(Qi,(32,---,(5n)- G*^ Gf' ... G 




iT 


where the Qi matrices are elements of GL{d,M.) and 
there is a pair (fc, 1) for which Qk ^ Qi. 

Now 


G **—I/O** _ /O* 

k — ^k 


^Qf^QiGt ^ Gr'Gr. 
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For the collection {G*}igv of matrices in GL{d,R), let 


Lemma 6 states that connectivity is a necessary prop¬ 
erty to determine a unique (up to transformation from 
the left) collection {G*} satisfying (2). As it turns out, a 
stronger type of connectivity - quasi-strong connectiv¬ 
ity - is useful in order to develop linear algebraic meth¬ 
ods for solving our synchronization problem. The first 
method we present is based on the so called Z-matrix. 

3 Linear invertible transformations 

3.1 The Z-matrix 

In this section a certain matrix is defined - referred to 
as Z. It is used as a building block in a matrix iJ, corre¬ 
sponding to the Hessian of a convex quadratic function, 
see Section 3.4. After its definition, its properties are in¬ 
vestigated. Amongst other things, it is shown that if the 
Gij transformations are orthogonal, i.e., GfjGij = I, the 
matrix (—Z) is (critically) stable in the linear dynami¬ 
cal systems sense (cf. Lemma 14). This means that, for 
directed graphs, the matrix Z can be used in a linear 
distributed algorithm for synchronizing orthogonal ma¬ 
trices (Section 4.1). 

Define the matrix 


C/i({G*}.ev) = [g*-^ g;-^ ... Gtr^Y, 
C^2({G*}iev) = [g^ GI...GI - 

Lemma 8 For any (QSC) graph Q = {V,S), collec¬ 
tion {GL}(j^j)gf of matrices in GL{d,M.) - transitively 
consistent for 0 - and collection {G*}igv of matrices in 
GL((i, R) it holds that 

GT = G*-^G* for all {i,j)eS 

(if and) only if 

im{diag{Gl,G 2 , ■ ■ ■, G*) V) = ker{L 0 Id), (3) 

for any matrix V, where the columns thereof form a basis 
for ker{Z{0,{G*j}(^ij)^£)). In particular, if Q is QSC, 
(3) can he stated as 

zm(C/i({G:}.ev)) = ))• 

Proof: 

Only if: Suppose it holds that 

Gb = G*-^G) for all {i,j) G 8. 



Then 

Z{Q,{G*,}(i,,)^e) (4) 

= diag(Gr\Gr\...,Gri)(L®/)• 

diag(Gt,G*,...,G:). 


and the matrix 


Now, 


Z{Q, {Gb}(ij)gf:) 

= diag(A(e)l) ®Id- W{g, {Gb}p,,)e£). 

The symbol 0 denotes the Kronecker product. 


Z{g,{GT}i^,j)ee)V = 0^ 

{L 0 /)diag(Gt, G^,..., G*JV = 0 
im(diag (GI, G^,..., G:) y) = ker(L 0 Id). 


Remark 7 A more general way of constructing the W- 
matrix and the Z-matrix with positive weights is as fol¬ 
lows. Replace the Wij in the definition ofW with OijWij, 
and replace diag{A{g)l) 0 Id in the definition of Z with 

diag{ E aij , ^ • ■ • , ^ ^ *^nj) Z) Id' 

jeMi jeM2 

The Oij are positive for all i,j. Equivalent results to all 
the results obtained for the Z-matrix in this section can 
also be formulated for the alternative Z-matrix with pos¬ 
itive weights. The alternative Z-matrix can be used in a 
distributed algorithm, equivalent to the one that will be 
presented in Section 4.1.1. 


If: This part is only proven for the case when the graph 
g is QSC. 

Since {Gb}p is transitively consistent for g, there 
is {G**}igv of matrices in GL{d,'R.) such that 

■}(..,)6£) 

= diag(Gr-\ G**-\ ..., G"-i)(L 0 I)' 

diag(Gr,Gr,...,Gr). 

Thus, the null-space of Z{g, {Gb}(ij)g£) is given by 
ker(Z(^l,{Gb}(i,j)gf:)) = im(y), 
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where 

y = diag(Gr'**, G 2 - 1 **,..., G-i**)([l, Id). 


Now, suppose (3) holds. Then 

diag (G*i, G;, ...,Gl)V = ([1,1,..., 1]^ 0 Id)Q, 

where Q is some matrix in GL{d, R). This means that 

diag (G^Gr'”, G;G 2 '**,..., GIG-^**) =In®Q, 

which implies that {G**}igv and {G*}igv are equal up 
to transformation from the left. By using Lemma 6 we 
can conclude that 

GT = G*-^G* for all { 1 , 3 ) £ £■ 


Remark 9 In Lemma 8, the relation 

im{diag (G*, G 2 ,..., G* ) V) = ker{L 0 Id) 

holds if and only if for any matrix V 2 , where the columns 
thereof comprise a basis for ker{L ® Id), there is a matrix 
Q such that 

diag{GlG;,...,Gl)V = V2Q. 

Remark 10 In Lemma 8, ifQ is connected but not QSC, 
it can hold that is the adjacency matrix of a QSC 
graph Q' = (y',£'). Then it holds that 

*m([/i({Gr^}.ev) = ker{Z{g',{G*J}^,,,)^s). 


Lemma 8 is important as it provides a way of finding ma¬ 
trices {G*}igv fulfilling (2). In the following subsection, 
this lemma is used to provide a least squares method. 

3.2 A least sguares method 

Suppose the graph Q = (V, £) is QSC, and the collection 
{Gij}(ij)^£ of matrices in GL{d,R) are not transitively 
consistent for Q, but close to being transitively consistent 
(closeness is in the sense of some matrix norm in 
Then, motivated by Lemma 8, the collection {G*}igv of 
matrices in GL{d,M.) such that (2) holds can be found 
by using the following approach. 


Algorithm 1 

(1) Solve the problem 

niin||Z(^l,{Gy}(,,j-)gf:)C|||, 

where V G R^dxd^ yTy ^ p,y 

means of the Singular Value Decomposition of 
Z{g, {Gij}(j j)g£). Let Vi be the optimal solution. 

(2) Identify the G* in the collection {G*}igv by 

Vi^ = [Gr^Gr^,...,Gr^]. 

The algorithm is motivated by Lemma 8. Note that the 
method is applicable if and only if the graph Q is QSC, 
(Lemma 8). In the special case when the transforma¬ 
tions are known to be Euclidean (or belong to some other 
desirable subset of GL(d, R)), the collection {G**} can 
be obtained by projecting the G* onto the set of Eu¬ 
clidean transformations (or any other desirable subset 
of GL(d, R)). 

If {Gij}(^i jy£ is close to being transitively consistent, 
the d X d block matrices in Vi are invertible and can 
be identified with the G^”^. This is guaranteed by the 
following lemma [23]. 

Lemma 11 In this lemma Z or Z{g,{Gij}(^i jy£) ks 

fixed, whereas the matrix Z is regarded as a variable in 
W^dxnd^ Let 

Si = {U e : U'^U = I}, 

S 2 (Z) = arg min tiace(U^Z^ZU). 

For e > 0, there is (1(e) >0 such that if 

\\z-zy<s, 

it holds that for all U G £ 2 ( 2 ), 

\\U\\s2{z) < e, 

where 

3.3 Further results 

Loops in the graph Q are essential for the performance 
of Algorithm 1 - if the graph is QSC and has no loops, 
improvement is not possible, see the following lemma. 

Lemma 12 If the QSC graph G is a spanning tree (con¬ 
taining a center), any collection {Gij}(^ijy£ of matrices 
in GL{d,M.) is transitively consistent for G. 
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Proof: 

Z{G, 

/ * ^ ... * * 
0 /*...* * 


0 0 0 .../* 

0 0 0 ... 0 0 

where all but one of the * at each (block) row is nonzero 
and an invertible matrix. Due to this structure, there 
is a collection {G*}i^v of matrices in GL{d,M.) such 
that (4) holds, which in turn means that Gij = G*~^Gj 
for all {i,j) G S. Now since G is QSC, this means that 
Gy =G*-^G'*forall(i,j)GVx V. ■ 

Due to Lemma 12, if G is QSC and a spanning tree and 
if corresponds to “disturbed” versions of 

{G*j}i^ij'j^£, the solution to Algorithm 1 will only pro¬ 
vide the once again. 

Lemma 11 provides us with the positive result that the 
solution to Algorithm 1 depends continuously on the Gy 
transformations. A somewhat negative result is provided 
by Lemma 13 below. Unfortunately it is not true that 
(3) implies transitive consistency. 

Lemma 13 Let G = (V, £) be any QSC graph satisfying 
that at least one element in the vector A{G)[l, 1, • ■ •, 1]^ 
is greater or equal to 2. Let {G*j}{ij)(z£ he a eolleetion of 
matriees in GL{d,M.), transitively eonsistent for G ■ Let 
{G* }igv be a eolleetion of matrices in GL{d, R) for which 
it holds that 

GL = G*-^G* for all (z,j) G f. 

Now, for any e > 0, there is a eolleetion {Gy}(jy)g£ of 
matriees in GL{d,M.) that are not transitively eonsistent 
for G sueh that 

E \\G^J-GL\\p<e, (5) 

and (3) holds for {Gij}(ij)^£ and a collection {Gijigv 
of matrices in GL{d, R.). 

Proof: 

Suppose the fcth element of the vector A{G)[1, 1, ■ ■ ■, 1]^ 
is larger or equal to 2. Then there is I, m such that I ^ k, 
m ^ k, GJyGJ^ G GL{d,R). For 0 < a < 1 let Gki = 
(1 + a)Gli and Gkm = (1 - <^)GX^. Furthermore, let 


Gy = G*j for all {i,j) ^ {{k,l),{k,m)}. It is easy to 
see that the left-hand side of (5) is less than or equal to 
adlGfcillF + IIGJ^IIf). Now we choose 

e 

and (5) is satisfied. By construction, all the Gy are ele¬ 
ments of GL{d, R). 

Let Gi = G* for all i. It holds that 

Z{G, {Gy}(iy)gf:) 

= dmg{Gf\Gf\ ..., Gf^){{L + Q)® /)• 
diag(Gi,G 2 ,.. ■ ,G„), 

where Q = [Qij], Qki = a, Qkm = -ct and Qy = 0 
for all {i,j) ^ {(fc,0; ik,m)}. Since ker((L + Q) ^ I) A 
ker(L 0 /), (3) holds for the Gi. According to Lemma 8, 
if the Gy are transitively consistent and G is QSC, (3) 
is a condition to guarantee (2). But (2) is not fulfilled 
since GkGuGj^ = (1 -I- a)I. Thus, the Gy are not 
transitively consistent. I 

After the introduction of Lemma 13, one might be lead to 
believe that Algorithm 1 does not work well in practice. 
However, as will be seen in Section 5, this is definitely 
not the case. 

Now, to recap: Transitive consistency is equivalent to 
(2). Lemma 8 states that when G is QSC and transitive 
consistency holds, (2) and (3) are equivalent. However, 
Lemma 13 states that (3) is not equivalent to transitive 
consistency for QSC graphs. 

Now we show a stability property oi —Z. If the Gij are 
transitively consistent, it is easy to see (from (4)) that 
— Z(G, {Gy}(iy)gf:) is critically stable, see definition in 
Lemma 14 below. However, the following result shows 
that if the Gy transformations are elements in 0{d), i.e., 
GfjGij = I for all i,j, the matrix -Z{G,{Gij}(^ij)(z£) 
is critically stable regardless if transitive consistency is 
fulfilled or not. 

Lemma 14 For any graph G = iV,£) and collection 
{Gy}(iy)g£ where Gij G 0{d) for all (i,j) G £, the 
matrix —Z{G,{Gij}(^i jf^£) is eritically stable, i.e., for 
any e > 0, there is (5(e) such that for a;(0) = xq G R”‘^, 
||a;(0)|| < S it holds that 

lk(Q|| < e, 

when 

x(t) = -Z{G,{Gij}(^ij)^£)x{t). 

Furthermore, if there are eigenvalues exactly on the imag¬ 
inary axis, those eigenvalues are equal to zero. 




Proof: 

Let 

x(t) = —Zx{t), x{t) G (6) 

where a;(^ is the initial state. We can write x{t) as 
x{t) = [x'{{t),x'^{t),... ,x'^{t)Y', where Xi{t) G for 
all i. Define the function 

V{x) = max(xfxi). 

i 

If there is some eigenvalue of Z with negative real part 
or if there is a Jordan block of dimension larger than one 
corresponding to an eigenvalue on the imaginary axis, 
there is xq such that for the state x{t) with initial state 
xq, V{x{t)) —>■ oo as t —>■ oo. We want to show that this 
is not possible. Let us first define the set 

21max(i) = {* : v{x(t)) = xf(t)xi{t)}. 


Now, 


D+{V{x{t))) = max —xf{t)x^{t) 
lelmaxl*) dt 


(7) 


= . max xj{t) V {G^JXJ{t) - x^{t) 


< 0 , 


where is the upper Dini-derivative. A proof of the 
first equality (7) can be found in [29] using the results in 
[30] and [31]. The result appears frequently in the litera¬ 
ture [32,33]. Now we can use the Comparison Lemma [34] 
to show that V{x{t)) is decreasing independently of the 
choice of xq. The inequality in (7) is a consequence of 
the fact that the Gij are orthogonal matrices. 

Now we show that there are no non-zero eigenvalues on 
the imaginary axis. Suppose there are non-zero eigen¬ 
values on the imaginary axis, then there must be a 
nontrivial periodic solution x{f) = [xf, xj,..., x'^]'^ 
to (6), i.e., x{t) is periodic and x{ti) ^ x{t 2 ) for some 
ti ^ t 2 - It can be shown that D'^{V{x{t))) = 0 for all t 
and it can also be shown that a necessary condition for 
this to hold is that Xi{t) = Xi{t) for all t and Gij = I 
for all The procedure to show the latter is a bit 

intricate and is based on an induction argument hing¬ 
ing on the fact that Q is QSC. Now, if the Gij ^ /, 
the necessary condition is not fulfilled, hence we have 
a contradiction. In the case when the Gij = / it holds 
that Z{Q) = L{Q) 0 Id and the latter matrix does not 
have any non-zero eigenvalues on the imaginary axis. I 


3-4 Optimization problems and the H-matrix 

In this subsection a matrix H is defined as the Hessian of 
a quadratic convex function. In the previous subsection 


the approach was to define a set of linear constraints, 
which are fulfilled for transitively consistent transfor¬ 
mations, and then use these constraints to formulate a 
least squares optimization problem. In this section the 
approach is different. Optimization problems are formu¬ 
lated directly, without taking a detour via algebraic con¬ 
straints. An assumption throughout this section is that 
Q is connected. 

Given the graph Q = (V,£l) and the collection 

{Gij}(ij)^£ of matrices in GL(d, R), we formulate three 
optimization problems, where the the first, (PI), cor¬ 
responds to the exact problem we want to solve. The 
objective function is non-convex and the constraint set 
is non-compact. The second problem (P2) is a restric¬ 
tion of the first problem having a compact constraint 
set (with a non-convex objective function). In contrast, 
the third problem (P3) has a quadratic convex objective 
function of the G~^ as well as a compact constraint set. 


(PI) 


min X; \\\Grj-GGGj\\l, 
s.t. Gi G GL(d,R). 


(P2) 


min X; ^\\G^j-GGGJ\\%, 
{(^iUev {i,j)£e 

s.t. Gi G 0{d). 


(P3) 


min E 

s.t. Gi({Gi}igv)^Gi({Gi}igv) = <5 0. 


Define the two functions 

f{{G-%^v)= E l\\GvG-^-G-^\\l, 

5({GJ.ev)= E IWGi, - G-^G,\\j,. (8) 


The matrix Q is symmetric and positive definite. We 
implicitly assume that g and / are parameterized by 
{Gij}(ij)^£. There is a similar problem to (P3), de¬ 
fined by left-multiplication by the Gi instead of right- 
multiplication by the G~^: 

[ min E k\\<^^Gij - GjWj,, 

(P4) ^ {Giliev (ij)gf: 

[ s.t. G 2 ({Gi}igv)G 2 ({Gi}igv)^ = Q 0, 

The two problems are equivalent. We choose to study 
(P3) instead of (P4) in order to more easily see the con¬ 
nection between the Hessian (the J7-matrix) in the prob¬ 
lem and the matrix Z (cf. Section 3.5). 
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The problem (P2) and variants thereof has received at¬ 
tention lately [1]. Exact solutions do not exist in gen¬ 
eral and local gradient descent methods are used. One 
of the more important contributions of this work is that 
we provide a lower bound for the global solution of this 
problem as well as a method for which the bound is al¬ 
most tight in numerical experiments. 

3.5 Problem (PS) and its connection to problem (PI) - 
definition of the PI-matrix 

Let X = C/idGijigv)- By a slight abuse of notation, let 
f{X) = f{{G-%^v)= E \\\GlG-^ - GfWl- 

(hj)e£ 

Now 

where 

= Z{g,{Gij}(^ij')^£) + Z2{G,{Gij}(ij)^£); 

and 

Z2{G, {Gij}(ij)^s) 

= diag{W{G,{Gij}(^i j-j^£)W{G, {Gij}(^ij)^£)^) 

-W{G, 

Gij = Gjj for all i,j and G = (V,f) is the graph 
constructed reversing the direction of the edges in G- 
The operator diag(-) is here understood in the block- 
matrix sense, i.e., for a matrix B G diag(i3) = 

{In ® IrflJ) © B, where 0 denotes element-wise multi¬ 
plication, In is the n-dimensional identity matrix and Id 
is the d-dimensional vector containing ones. 

Remark 15 A more general formulation of the objective 
functions f and g with positive weights is 

/({G-d^ev)= E lai,\\G.,G-^ - G-Yf, 

d,3)e£ 

giiGijiev) = ^ 2^i3\\Gij - GGGj\\%. 

(i.pes 

The Oij are positive for all i,j- This way of defining the 
objective functions lead to a slight modification of the H- 
matrix. Eguivalent results to all the results obtained in 
this section for the H-matrix can also he formulated for 
this alternative definition of the H-matrix with weights. 
The alternative H-matrix can also be used in an equiv¬ 
alent distributed algorithm to the one presented in Sec¬ 
tion j.1.2. 


Lemma 16 In the special case when all the Gij are ele¬ 
ments ofO{d), i.e., orthogonal matrices, 

Z2{G,{Gij},^ij'j^£) = Z{G,{Gij}(ij)^£). 

Furthermore, if the graph G is also symmetric, 

Z2{G,{Gij}(^ij^^£) = Z{G,{Gij}(^ij)^£)'^. 

Lemma 17 For any connected graph G = {V,S), and 
collection of matrices in GL{d,M.), the col¬ 

lection {G*^}(ij)g£ is transitively consistent for G if 
and only if there is a collection {G*}igv of matrices in 
GL((i, R) such that 

im{Ui{{G*}i(.v)) C fcer(iJ(0,{GL}(ij)g£:)). 


The collection {G*}igv satisfies (2). 

Proof: Suppose {Gy}(j is transitively consistent, 
then, according to Lemma 6, there is {G*}igv such that 
(2) holds, which in turn can be used to show that 

/({Gr'}.ev) 

= [/l({Ga^6v)^i^(e,{G*,}p,,)e£)C/l({Gd.ev) 

= 0. (9) 

On the other hand, if {G*j}(^ij)^£ is not transitively 
consistent, there is no {G*}igv such that (2) holds. It 
can now be shown that (9) does not hold for any collec¬ 
tion {G*}igv of matrices in GL{d,M.). ■ 


Lemma 18 For any connected graph G = {V,S) and 
collection {G*j}pi j'^(z£ of matrices inGL{d,M.) - transi¬ 
tively consistent for G - it holds that 

dim{ker{H{G, {GL}(ij)gf:))) = d. 


Proof: Due to Lemma 17, we know that 

dim(ker(iJ(e, {GP}^ij)e£))) > d. (10) 

Thus, we need to show that the inequality in (10) can¬ 
not be strict. Since |G*,|c, is transitively consistent, 

there is {G-}.evMfiili.)(d 

Suppose the inequality is strict for {G*j}(^ij)^£. We know 
there is {G*}igv where G* G GL{d, R) for all i, such that 

im{U,{{G*}^ev)) C ker(i7(e,{GL}(,,,)e£)). 
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Now there must be a vector y = [yf, yj,..., y'^]'^ G 
where the yi are in such that 

yeker(i7(c;,{Gb}(,,,.)e£)), 

y 0, and y^C/i({G*}igv) = 0- There must be k and 
I such that the ^th element of yk is nonzero. The set 
of transformations {G^“^G*}igv satisfy (2) (Lemma 6) 
and /({(Gr^G*)-i}iev) = 0. Now, let 

X = [5i,x2, ...,Xd] = Gi({G*-'G*},ev), 

and 

Y = [xi,X2,. ■. ,xi-i,y,xi+i,xd] = Ui{{Gl~^G*}i(zv). 

We know that H{0,{Gij}(ij)^£)Y = 0. For all i, let 
Yi be the zth d x_d block matrix in Y. We know by 
construction that Yk G GL{d,R). Now, for any j G A4 
it holds that 

\\Gl^Y,-Yk\\F = 0, 

which implies that Yj G GL(d,R). Also, for any i such 
that k G JVi, it holds that 

IIG:kYk - Yill^ = 0, 

which implies that Yi G GL(d,R). Now, due to the fact 
that is connected, an induction argument can be used 
to show that all the Yi are elements in GL{d, K.). 

The collection {Yi}i^v satisfies 

Gb = for all (z, j) G £. 

It_is easy to see that the two collections {Y~^}i^v and 
{G*} igv are not equal up to transformation from the 
left. But, since the graph is connected, the two must be 
equal up to transformation from the left (Lemma 6). 
This is a contradiction. Hence it is a false assumption 
that the inequality in (10) is strict. ■ 

We summarize the results of Lemma 17 and Lemma 18 
in the following proposition. 

Proposition 19 The collection {Gij}(^ij)^£ is transi¬ 
tively consistent for the connected graph Q if and only 
if there is a collection {G*}igv of matrices in GL{d,R) 
such that 

z7?z(C/i({G*}i6v)) = fcer(iJ(0,{GT}(ij)gf:)). 

The G* satisfy (2). 

Proof: Direct application of Lemma 17 and Lemma 18. 


The following proposition provides a similar, but some¬ 
what stronger result. 

Proposition 20 The collection {GT} (^ij)^£ of matrices 
in GL{d,R) is transitively consistent for the connected 
graph Q if and only if 

dim{ker{H{g, {GT}(,j)gf:))) = d. 

Proof: 

If: Let Y = [z/i, z/ 2 , • ■ •, Vnd]^ G be any full rank 

matrix such that 

H{g,{G:,}i.,,)e£)Y = o. 

All thejji G Let l) be the zth dxd block matrix in Y. 
Since Y is full rank, there is a (finite) sequence 
such that [z/Zi ,yi^,... yif\ G GL{d, R). 

Now, for fc G V we know that for any j G Afk it holds that 

\\Gl^Y,-Yk\\F = 0, 

which implies that im(l^T) = im(l)?’). Also, for any i 
such that k G Afi, it holds that 

which implies that im(FT) = im(i)J). Now, due to the 
fact that Q is connected, an induction argument can be 
used to show that im(l^T) = ini(FT) for all z, j. But then 

im([z/*i, z/* 2 , • ■ • yid]) C im(y/) for all j, 

which together with the fact that [z/Zi, z/Za) • ■ • G 
GL{d,R) is full rank, can be used to show that 
Yi G GL{d, R) for all i. Thus, 

im(Gi({y-i}zev)) = ker(iJ(e, {GL}(,,,)e£)), 

and the desired result follows from Proposition 19 where 
the G* are replaced by the Y~^. 

Only if: Direct application of Lemma 18. ■ 

Lemma 21 The optimal solution to (PS) is 
X* = VP, 

where P G andV G The matrix V is given 

by the solution to the problem 

imn trace(lF^i7(1/, {Gy }(z j))IT), 

W G = I. 
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x{t) = -i?(^?,{Gy}(ij)gf:)x(i). 


and P is given by the solution to the problem 


when 


(P6) 


min trace{P^{V^H{g, {G,J}i^,J)es)V)P), 
P G ^ Q 


Proof: In the new notation, problem (P3) is written as 

( mintrace(P^(IP'^iJ(0,{Gjj}(j_j)g£)IP)P), 
w,p 

[ p(Z pTp (Z ]gndxd^ W^W = I. 

In the following derivations it is assumed that P and 
W belong to the constraint sets defined in the problem 
above. 

mintTace{P'^{W'^ H{g , {Gij}(^i j)(zs)W)P) 

w,p 

= mintTace{P'^{Q'^{W)D{W)Q{W))P) 

W,P 

= mintrace(P^P(IP)P) 
w,p 

= nnntrace(P'^(I/^P(0, {G^j}pp(,£)V)P), 

where (W)D{W)Q{W) is the spectral factorization 
of 

W^H{g,{G,,}p,)^e)W. 


Proposition 22 For any Qi >- 0 and Q 2 P 0 let 
{G*}igv <^nd {G^ligv be the transformations obtained 
from the optimal solutions of problem (PS) with Q equal 
to Qi and Q equal to Q 2 respectively. It holds that 


Proof: 

The matrix H{g,{G*j}(^ij')^£) is the Hessian matrix 
and hence positive semi-definite. I 


3.6 A least squares method 

Proposition 22 is important, it states that we can with¬ 
out loss of generality assume that Q = I, since the choice 
of Q does not affect the value of g, i.e., the cost function 
we want to minimize. The value of / changes with Q, 
but this is of less importance. Motivated by these results 
we introduce a least squares method along the lines of 
Algorithm 1. 

Algorithm 2 

(1) Let Vi be the optimal solution to problem (P5). 

(2) Identify the G* in the collection {G*}igv by 

ifi^ = [Gr^Gr^,...,Gr^]. 

The algorithm is motivated by Proposition 19 and 
Proposition 22. If the collection {Gij}(ij-)^£ is close 
enough to be transitively consistent, step 2) can be 
executed, i.e., each d x d sub-block of the matrix Vi is 
invertible. The result that guarantees this is analogous 
to the statement in Lemma 11. 

3.7 A Gauss-Newton method 


g({G*}^6v) = ff({GrW), 

i.e., the value of g is independent of Q. 

Proof: According to Lemma 21 the transformations are 
equal up to transformation from the left. I 


Remark 23 It is implicitly assumed in Proposition 22 
that the G* and the G** are in GL{d, R). This is guaran¬ 
teed if the Gij are sufficiently close to be transitively con¬ 
sistent. The result to guarantee this is omitted but anal¬ 
ogous to the statement in Lemma 11 for the Z-matrix. 

Lemma 24 For any graph g = {V,S) and collection 
{Gij}(i,j)^£ where Gij G GL{d,M.) for all {i,j) G S, the 
matrix —II{g,{G*j}(^ij')^£) is critically stable, i.e., for 
any e > 0, there is d{e) such that for any a;(0) = xq G 
||a;(0)|| < 6 it holds that 

IkWII < e, 


In this section a Gauss-Newton method is presented. The 
solution obtained in Algorithm 2 is used as the initial¬ 
ization for the algorithm. 

The Frechet derivatives of the identity map and the in¬ 
verse map at the point Gi in the direction Ei are given 

by 

Lid(Gi, Ei) = Ei, 

Linv{Gi, Ei) = —Gj ^EiG^ 

respectively. Higham [35] provides a good introduction 
to Frechet derivatives for matrix functions. Let {Pi}jgv 
be a collection of matrices in R"^”. It holds that 

g{{Gi -\- Pijigv) 

= E l\\G^j-GpG,-GpLi£i{G„Ep 
— T\nv{Gi, Ei)Gj -f o(|| [Ei, Ej] 11^)111'. 
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Let 


When i = j, Hu is defined by 


= E \\\G.3-G-^G,-G-^Ua{G„E,) 

(ij)e£ 

~ Lin^(Gi, Ei)Gj\\% 

= E l\\G^j-G-^G,-G-^E, 

(ij)e£ 

+ G-^E,G-^G,\\l. 

Consider the following problem 

(P7) I 5({G*}iev,{^^i}iGv)- 

I {Ei}i€V 

Problem (P7) is solved in each Gauss-Newton step of the 
method we present below (Algorithm 3). Its solution is 
given by the collection {Ll*}iev obtained by 

vec(U2i{E*}iev)) = X, ( 11 ) 

where x is obtained by the solution of 

-ffGN({Gi}igv, {Gij}(ij)^£)x 

= — CGN({Gi}igv, 5 {G'ij}(j^j)g£:); (12) 

vec(-) is the vectorization operator, i.e., it returns a vec¬ 
tor with the stacked columns (in consecutive order) of its 
matrix-argument. The matrix iJcN G and the 

vector CGN G are defined as follows (for simplicity 
we have omitted the explicit dependence of {Gijigv and 
{Gij}{i,j)e£)' 

— \Eij ], 

where Hij £ for all i,j. When i ^ j, Hij is 

defined by 


= E ((G-'G,GJG-^) 0 (G-^G-i))+ 

E Ici®{GfGf). 

{jiiGAbl 

Now, CGN = [cf, C 2 ,..., c^]^, where Ci G for all i. 
The Ci are defined by 

Ci = E ((G.-'G,) 0 G-^)vec(Gi, - G-^G,)- 

jeHi 

Y, Id® G-^vec(G,i - G-^G,). 

{j-.i&Mj} 

Algorithm 3 

(1) Run Algorithm 2 and let {G*}igv bet the collection 
of matrices obtained in step (2) of that algorithm. 

(2) Let 9 G* = 0 for i = 1, 2,..., n. 

(3) repeat: 

(a) Gi ^ Gi + E* for all i, 

(b) Update the E* by (11) and (12), i.e., 

vec(G2({Ll*}i6v)) = X, 

where x is the solution to (12). 

The sloping criteria in step (3) of Algorithm 5 could 
be that the improvement of the cost function is smaller 
than a certain threshold for two consecutive iterations, 
or it could be that a certain number of iterations have 
been executed etc. It should be noted that TLgn is both 
positive definite and sparse. In order to solve (12) one can 
use for example the Conjugate Gradient method [36,37]. 

3.8 Affine and Euclidean transformations 


H^3 = 


{0 


f(j ^M), 
[{i^Afj), 

{-((G-iG,)®(G-^G-i)) 


f(j gM), 

{-((GfG-^)®(G-^G-^)) 

it] 

f(j ^M), 

l(* G AGO, 

r-((G-iG,)®(G-^G-i)) 


f(j GAG), 

l-((GfG-^)®(G-^G-i)) 

l(* G AG). 


In this subsection we consider affine and Euclidean trans¬ 
formations. These transformations are linear when ho¬ 
mogenous coordinates are used. To be more precise, an 
element in Aff((i, R) is a matrix 


G = 


Q t 
0 1 


where Q G GL{d, R), t G and 1 is a scalar. Its inverse 
is given by 


G-^ 


Q-^ -Q-H 

0 1 
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Euclidean transformations, E{d), is a special case of 
affine transformations where the matrix Q G 0{d). 

For any connected graph Q = (V, S) (due to Lemma 6), if 
and only if the collection {G'*j}(ij)g£: is transitively con¬ 
sistent and contains only affine transformations, there is 
a unique (up to transformation from the left by affine 
transformations) collection {Gijigv of affine transfor¬ 
mations such that 


ij — ^j' 


Each Gi is given by 


Qi ti 

0 1 


and each G*j is given by 


^ij — 


Qi ^Qj Qi \tj 

0 1 



Let 


where Hij G 
defined by 


EaU — \.^ij ] 5 

-i)x(d-i) £qj. j When i ^ j, Hij is 




{0 

‘7 


{-Q7^Q7^ 

•7 

f(J eM), 


•7 

l(iGA7j), 

l-Q7^Q7' 

•7 

f(J gM), 

1- Q7^Q7' 

l(* G A)-). 


When i = j, Hu is defined by 


j&Ht {j-ieAfj} 


Now, let {G*j}(^ij')^£ be a collection of matrices in 
Aff((i, M) that are not necessarily transitively consistent. 
It holds that (by a slight abuse of notation) 

ff({GJigv) = E ~ 

{id)e£ 

= g{{Qi}i(^v) 

+ E \\\id-Q^\Q-U)\\l, (13) 


where tij is the translational part of the transformation 
Gij. We see that there is a special structure of (13), 
where the cost function consists of two parts. The first 
part is only a function of the Qi , whereas the second part 
is a function of both rotations and translations. 


The matrix Hag and the vector cag depend on {Gjjigv 
and {Gy . 

The solutions to the problem (P8) is given by the ele¬ 
ments in the set 

{{U}i(zv : , tlf = -cag}- 


The Gauss-Newton method developed in Section 3.7, 
i.e., Algorithm 3, can be adapted to the case of affine 
transformations. Now we require that 


Ei Q B = Ei for all i, where B = 


IdlJ Id 
0 0 


Define the following optimization problem 


(P8) 


[ min Y} 

y {ti}iev {iu)e£ 



Qi -1*)Ill- 


Let 


CAff \c.-^ 7 ^2 5 ■ * ■ 5 \ 1 

where 

Ci = Qj tij — Q^ tji- 

j&Hi {j-ieJdj} 


We remind the reader that 0 denotes element-wise mul¬ 
tiplication. In each iteration (in the modified step (3) of 
Algorithm 3) the collection {E*}igv is obtained by 


vec([/2({E*}igv)) = X, (14) 

where x = Xv, v is obtained by the solution to 

X^HgnXv = -A^cgn, (15) 


and A G ig defined below. 


A = J„ 0 (/d 0 B), 
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where 

B = G ^dxid-i)^ 

0 

Now we present the following algorithm for affine trans¬ 
formations. 

Algorithm 4 

(1) Run Algorithm 2 for the collection {Qij}(ij)(£S 

let {Qi}i£v be the matrices obtained in step (2) of 
that algorithm. 

(2) Solve problem (P8) for the ti using the Qi from (1) 
and let 

Gi = * for all i. 

_0 1 _ 

(3) Let 9 A* = 0 for f = 1, 2,..., n. 

(4) while a stoping criteria has not been met: 

(a) Gi ^ Gi + E* for all i, 

(b) Update the E* by (14) and (15), i.e., 

vec([/2({Ll*}*6v)) = X, 

where x is the solution to (14). 

Remark 25 There are many variations of Algorithm 4 
that can be employed. The most simple one is to omit 
steps (3) and (4). Another one is to run the Gauss- 
Newton method (Algorithm 3) for the Qi matrices after 
step (2). The expression in (13) can also be changed to 
include weights. For example, if the orthogonal matrices 
are closer to be transitively consistent than the transla¬ 
tions, the first part of the expression, i.e., g{{Qi}iev)> 
could be weighted with a positive weight larger than 1. 

Remark 26 After a slight modification. Algorithm 4 can 
be used for Euclidean transformations instead of affine 
ones. In this case Algorithm 5 (see Section 4) is used 
in (1) to generate the Qi transformations instead of Al¬ 
gorithm 2. Numerical simulations (see Section 5) show 
that this is a good method in comparison to Algorithm 1 
or Algorithm 2 (where the matrices are finally projected 
onto the set of Euclidean transformations E(d)). 

4 Orthogonal matrices 

In this section problem (P2) is studied. For orthogonal 
matrices the objective functions / and g are equivalent. 
The Gauss-Newton method (Algorithm 3) is hence not 
necessary. Furthermore, the orthogonal matrices is an 


important class of matrices, not the least in dimension 
d = 3. 

We begin by formulating the following result. 

Proposition 27 For the connected graph G = (V,£l), let 
{Gij}(ij)^£ be a collection of matrices in GL{d,M.). Let 
{Gijigv be a collection of matrices obtained from Algo¬ 
rithm 2. Let {G*}igv be a collection of matrices solving 
the optimization problem (P2). It holds that 

/({v^G-H.ev)<5({G*} iev)- 


Proof: 

It is easy to verify that for orthogonal matrices, (P3) is a 
relaxation of (P2) when Q = nl. Now, (Proposition 22) 
the solution to (P3) with Q = nl is provided by the ma¬ 
trices obtained by Algorithm 2 after scaling by ■ 

Let us now extend Algorithm 1 (Algorithm 2) in the 
following way. 

Algorithm 5 

(1) Same as in Algorithm 1 (same as in Algorithm 2). 

(2) Same as in Algorithm 1 (same as in Algorithm 2). 

(3) Let G** be the projection of G* onto 0{d), i.e., 

G** = QV^, 

where QDV^ is the SVD of G*. Let 

^ij • 

The collection is the final transitively 

consistent collection. 

Proposition 27 can now be used to provide performance 
guarantees. An upper bound on the closeness to opti¬ 
mality is given by 

5({Gr}iev))-/(v^Gr'W), (16) 

where the G* are obtained from Algorithm 2 and the 
G** are obtained from Algorithm 5 - assuming the first 
two steps are the same as in Algorithm 2. 

If the Gij are also elements in 0{d), the difference in (16) 
is almost tight. For example, in the case when d = 3, 
n = 100, and the Gy are generated from Gi-matrices 
and -matrices matrices by Gy = Gf^GjRij {Rij is 
an orthogonal matrix with geodesic distance to / less or 
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equal to 7r/4. It is generated by drawing a skew symmet¬ 
ric matrix from the uniform distribution over the closed 
ball with radius 7r/4 and then taking the matrix expo¬ 
nential of that matrix). Let 

M{Grw,{Gr'}. 6 v) 

_ g({Gt*}igv)) —/({y^G* ^}igv) 

fiiVnG* ^}iev) 

For 1000 experiments we observe that 

M{Gr}igv,{Gri}.6v)<6*10-4. 

This means that the solution obtained by Algorithm 5 
is closer than 0.06% to the global optimum of problem 
(P2). The graphs in these experiments were QSC and 
the adjacency matrices contained 100 zero entries. 

4-1 Distributed algorithms 

In this subsection we show that Algorithm 5 can be 
implemented in a distributed way. Besides the graph 
G = which describes what transformations are 

available, another graph is used. It is 

always assumed S C The graph is referred 

to as the communication graph. The assumptions on the 
commnnication graph differ between the two pre¬ 
sented algorithms. 

4-1.1 Orthogonal matrices and QSC communication 
graph 

Here it is assnmed that all transformations are orthog¬ 
onal matrices, i.e., elements in 0{d). That is, the Gij 
matrices as well as the G*, matrices and the G* matrices 
are assumed to be elements in 0{d). 

The algorithm will now be presented, after which an ex¬ 
planation and justification is provided. In this algorithm 
it is assumed that Q = is QSC. The notation A/) is 
used to denote Mi(Q) = 


Algorithm 6 
Let 

A(0) = [An0),Aj(0),...,Aj(0)f, 

where, for all i, the elements of the matrix Ai(0) G 
are drawn from W(—0.5, 0.5), i.e., the uniform distribu¬ 
tion with the open interval (—0.5, 0.5) as support. Let 
Xi{t) for < G N be defined by the following distributed 


algorithm: 

Xi{t + 1) = Ai(t) 4- e ^ {GijXj{t) - Xi{t)), 
jeM 

X2{t + 1) = X2{t) + e ^ {G2jXj{t) - X2{t)), 

J 6 AA 2 


A„(t + 1) = A„(t) + e ^ {Gn,X,{t) - A„(t)), 

where e > 0.1 In compact notation this is written as 

X{t + 1) = Xit) - eZ{g, {G.,}(,,,)e£))A(t). (18) 

For a sufficiently large t, let Gf* be the projection of 
Xi{t) onto 0{d), and let G*j = Gf*G* for all i,j. It 
should be noted that if the spectral radius is not known, 
in practice it is enough to choose e to something small. 

Analysis of the algorithm 

In this section the theoretical analysis of the algorithm 
is provided. The first thing we need to guarantee is that 
the matrix 

I - ^Z{g, {Gij}(ij)gf:), 

appearing in the right-hand side of the discrete time lin¬ 
ear system (18), is critically stable in the linear dynami¬ 
cal systems sense. This means that all eigenvalues must 
be smaller than or equal to 1 in absolute value and any 
Jordan-block corresponding to an eigenvalue whose ab¬ 
solute value is 1 must be one-dimensional [38]. 

Lemma 28 Ifg{V,S) is QSC and e > 0 small enough 
it holds that 

I - eZ{g, {Gij}(ij)g£:) 
is critically stable. 

Proof: According to Lemma 14 it holds that Z is crit¬ 
ically stable and has no non-zero eigenvalues on the 
imaginary axis. This means that for e small enough 
the eigenvalues of eZ{g, {Gij}(ij)g£) are located in the 
closed unit disc centered at —1; the eigenvalues on the 
boundary are simple. I 

Remark 29 Numerical simulations seem to indicate 
that in practice one can choose 

p{Z{g,{G.AcA<^e))) ^ 

where p{Z{g, {G^}(jj)g£)) is the spectral radius. 
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Now we can deduce that if e is chosen small enough, X (t) 
converges (to something). It easy to verify (Lemma 8) 
that if the Gij were transitively consistent, X{t) would 
converge with exponential rate of convergence to 

X = [xi,X2,...,Xd] = [X'[,X2,---,Xnf, 

where Xi is the projection of the tth column of ^(0) 
onto ker(Z) and Xi G for all i. Since the Xi(0) 

are drawn from the distribution (W(—0.5,0.5))'^^'^, it is 
extremely unlikely that (probability zero) X has not full 
rank. If all the Xi are full rank matrices, 

Gij = XiXj for all i,j- 


Now, if the Gij are not transitively consistent, in general 
the Xi{t) converge to 0, which is not favorable. How¬ 
ever, if the Gij are close to being transitively consistent, 
since the eigenvalues of Z are continuous in the Gij , the 
d smallest eigenvalues are significantly smaller in mag¬ 
nitude than the other eigenvalues; also the d smallest 
singular values are significantly smaller than the other 
singular values. Up to rotation, the right-singular vec¬ 
tors corresponding to the d smallest singular values are 
continuous in the Gy, see Lemma 11. 

Let the right-singular vectors corresponding to the d 
smallest singular values comprise the columns of the ma¬ 
trix Y € The matrix Y is equal to V obtained in 

the first step of Algorithm 1 (up to transformation from 
the left). Now, as t —>■ oo, under the assumption that 
the Gij are close to the G*j, the columns of X(t) con¬ 
verge to im(F) much faster than X(t) converges to 0. 
Thus, for t large enough X(t) is approximately equal to 
Y up to transformation from the left. This convergence 
can be seen in Figure 9 for different choices of n, d, and 
magnitudes of noise. 

The last step of the algorithm is justified by Lemma 8. 

4G.2 Orthogonal matrices and symmetric connected 
communication graph 

In this section a general distributed algorithm is pre¬ 
sented, which works for Gy matrices in GL(fi,M), a di¬ 
rected connected graph Q, and a symmetric communica¬ 
tion graph We will make the assumption that 
is the union graph of G and G, i-e., C/™™ = (V, £TJf). The 
difference between Algorithm 6 and Algorithm 7 pre¬ 
sented here, is that the Z-matrix is used in the former, 
whereas the iL-matrix is used in the latter. Here, differ¬ 
ent from Section 4.1.1, it does not hold that A/i(f/“™) = 
■hfiiG) for all i. When we write A/) this is shorthand for 


Algorithm 7 
Let 

x(0) = [xn0),xj(0),...,xj(0)]^, 

where the elements of the matrix W(0) G are 

drawn from G(—0.5, 0.5). Let Xi{t) for t G N be defined 
by the following distributed algorithm: 

Xi{t + 1) = Xi(t) + e ^ QijXjit) - Vi,Xi{t)), 

jeM 

X2{t + 1 ) = X2{t) + e Y, {Q2jX,{t) - V2jX2{t)), 


Xnit + 1) = Xn{t) + e Y (QnjXjit) - U„,W„(t)), 

where 

('^’p(i7(l/,{Gy}(,y),^)))’ 
p{H{G, {Gy }(iy)g£:)) is the spectral radius and 

Qij = Gij + GJ^ , 

Vij =Id + GjiGji. 

In compact notation this is written as 

Xit + 1) = X{t) - eH{G, (Gy }(,y)e£))X(t). (19) 

For a sufficiently large t, let G*~^ = Xi{t) and let G^ = 
G*~^G*j for all i, j. It should be noted that if the spectral 
radius is not known, in practice it is enough to choose e 
to something small. 

Remark 30 In the definitions of Qij and Vij, in the case 
when (j, i) ^ G, the symbol Gij should be interpreted as 
the matrix in containing only zero-elements. 

Analysis of the algorithm 

The convergence analysis of this Algorithm is analogous 
simpler than that of Algorithm 6. The eigenvalues of the 
iL-matri are real since the matrix is symmetric. Instead 
of using Lemma 14, Lemma 24 is used instead. 

4-2 Gradient flow for orthogonal matrices 

Under the assumption that all the Gy are elements in 
0{d), we here provide a method, which will be used for 
comparison to our earlier methods. Results, along the 
lines of the ones presented in this section, can be found 
in [39,40,41]. 
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For all i, define the cost functions 

g, : (0(d))" ^ R+ 
by 

g,(Gi,G2,...,G„) = ^ Aij\\G,G,,Gj - I\\l. 

jeM 

The overall cost function 

g : (0(d))” ^ IR+ 

is equal to (/, i.e., 

n 

5({G*},ev) = ^3*(Gi,G2,...,G„). 


The (negative) gradient flow on (0(d))^ of ^ is given by 

G.(t) = - ^ A,iiG.{t)G.,G,itff (20) 

-(G.(t)G.,G,(t)^))G.(t) 

- ^ A,{{G,it)GlG.{tff 

-iG,it)GlG,{tf))G.it), 

for all i GV. 


Now we present an algorithm, which improves on Algo¬ 
rithm 5. However, as will be seen in Section 5, this im¬ 
provement is marginal. 

Algorithm 8 

(1) Run Algorithm 5 and let {G**} be the orthogonal 
matrices obtained in step (3) of the algorithm. 

(2) Solve (20) numerically (for example by using ode45 
in Matlab) for a sufficiently large time interval 
[0, T] with the G** as initial conditions. 

(3) Let {Gi(t)“^Gj(t)}(ij)gvxv be the collection of 
transitively consistent matrices. 

Remark 31 In step (3) of Algorithm 8, if the Gi{T) 
are not elements of 0(d) (due to errors from numerical 
integration), they need to be projected onto 0(d). 

5 Numerical verification 

In our experiments, we consider Algorithm 8 first. Sub¬ 
sequently, the centralised Z- and H-matrix methods are 


evaluated for different configurations. Eventually, the 
analogous distributed methods are used in our simula¬ 
tions. In order to compare the methods, an assumption 
throughout this section is that the graph Q - describing 
what transformations are available - is QSC. 

5.1 Generating graphs and transformations 

For each of the following experiments, the collection 
{G*}"gy is generated by drawing random matrices in 
0{d) [17]. From that, the (full) set of transitively consis¬ 
tent matrices {G*j = G*“^G*} ijev is created. The noisy 
set of pairwise transformations {Gijjij^v is generated 
by adding element-wise Gaussian noise with zero mean 
and standard deviation a to each G*j. After adding the 
element-wise Gaussian noise, the matrix is additionally 
projected onto 0{d). 

Furthermore, a quasi-strongly connected (QSG) graph 
with graph density p - not mix up with the spectral ra¬ 
dius of a matrix - is generated in the following manner. 
For generating a minimum QSG graph Q = {V,S), two 
lists are used. One list keeps track of the nodes that 
are already considered, and one list keeps track of 
the nodes that have not been considered. By a mini¬ 
mum QSC graph we mean a QSC graph that is a (span¬ 
ning) tree, i.e., one with exactly n — 1 edges. Initially, 
we set V = {1,..., n}, £ = {(i,i) : i G V}, = {r}, 

where r G V is a randomly selected node, and CP = 
{1,... ,n} — {r}. Then, the following procedure is re¬ 
peated n—1 times: pick random nodes i G CP and j G 
CP, add the edge (i,j) to £, and update £p and £P 
accordingly. After n—1 repetitions a (minimum) QSC 
graph has been generated. At this point we store the edge 
set £ and call it £^^'". Next, random edges are added 
to £ until the the density of the graph is larger than or 
equal to p, which is defined below. 

We remind the reader that A is the adjacency matrix of 
Q with elements A^. The graph density p{Q) is defined 

by 

The intuition behind the graph density is that it is the 
proportion of the number of present edges in Q with 
respect to a fully connected graph (having edges) 
excluding the edges in £^^^. With that, p = 0 denotes 
a minimum QSC graph, whereas p = 1 denotes a fully 
connected graph. Generating random QSC graphs with 
different values of the parameter p allows us to consider 
different degrees of missing transformations. 

Using the graph Q with density p, the collection 
{Gy}(jj)g£ is the one that is eventually used for the 
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evaluation. In the simulations, for each individual sub- 
hgure the simulations have been performed with 100 
random sets of orthogonal transformations (the transi¬ 
tively consistent ones and the synthetically generated 
noisy versions thereof) and QSC graphs. Shown in the 
sub-figures is the mean of all runs. 

5.2 Algorithm 8 - orthogonal matrices 

For d = 3, Figure 1 shows upper bounds on the gap 
between the optimal value and the value of the objec¬ 
tive function obtained by two methods - Algorithm 5, 
green curve, and Algorithm 8 , blue curve. In Algorithm 
8 , the initial states are given by the solution to Algo¬ 
rithm 5. The ODE in (20) is solved numerically in Mat- 
lab by ode45. For each number of coordinate systems n, 
100 simulations are conducted and averages are shown in 
Figure 1. In each simulation a set of transitively consis¬ 
tent orthogonal matrices {G'*j}(j_j)gvxV ^'^e generated 
from a set of orthogonal matrices {G*}igv according to 
the description in Section 5.1 below. The graph G used 
in each of the experiments is the complete graph. 

For a single numerical experiment. Figure 2 shows the 
improvement of h when Algorithm 8 is used. One can 
see that Algorithm 5 generates matrices that are close 
to a local optimum of problem (P2). 

It can be seen that only a marginal improvement can 
be made using the significantly more computationally 
expensive Algorithm 8 . Due to the heavy computational 
burden, in the following simulations we omit Algorithm 8 
and focus on the methods based on the Z- and H-matrix. 



Fig. 1. Upper bounds on the optimality gap, i.e., h, for the 
solution of Algorithm 5, green line, and Algorithm 8, blue 
line. The graph is complete. 



Fig. 2. Improvement of h when Algorithm 8 is used. In this 
case when n = 15, d = 3, and the graph is complete. 

5.3 Centralized methods for matrices in 0{d) 

In this set of experiments we compare the H-matrix 
method, the Z-matrix method and a (naive) reference- 
based method, where the latter serves as baseline for the 
comparison. 

5.3.1 The reference-based method 

For the reference-based method, a minimum QSC graph 
^min-QSC ^ £;mm-QSC g £■) ^ gjggg jg 

domly drawn as a subgraph of t/ = {V,S). For that, all 
centers (see Def. 2) of the graph Q are initially deter¬ 
mined by looking at the n x n distance matrix between 
all n nodes. From the set of centers, a node c is randomly 
selected. Since G is QSC, there is at least one such cen¬ 
ter. Let = {c} be the list of nodes that have al¬ 
ready been considered and initialise _ 0 ^ 

following procedure is repeated until |£:™™-QSC| _ 

randomly select a node r G CP , select a random 

node r' S {i : (i, r) G .?}, if there is such an r', add the 
edge (r',r) to ^““-QSC and add r' to . 

Per construction, the graph ^™in-QSC jg a spanning tree 
that contains a center. Thus, according to Lemma 12, 
the set {Gy }(i j)g£min-Qsc is transitively consistent for 
^min-QSC^ W.l.o.g., by Setting G* = / for the center c of 
^mm-QSC^ all (other) G* are (uniquely) determined as 

G* = G*G-.i for * p c, (z, j) G , ( 22 ) 
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To summarise, in the reference-based method a (ran¬ 
dom) rooted spanning tree graph is considered as sub- 




graph of Q, i.e., all but n — 1 relative transformations 
Gij (accounting for the transitive inconsistency) are dis¬ 
carded such that the remaining n — 1 relative transfor¬ 
mations are transitively consistent. 

In Fig. 3 the results of the experiments are shown. On 
the vertical axis, a normalised version of the function in 
(8), defined by 

5'({GJ.ev) = ^ E \\Gv-G-^Gj\\l (23) 

is used. Each sub-figure shows a different varying param¬ 
eter on the horizontal axis. The title of each sub-figure 
indicates the fixed parameters. 

It can be seen that in all cases the Z-matrix approach 
is nearly as good as the H-matrix approach when look¬ 
ing at orthogonal transformations. However, as antici¬ 
pated, the reference-based method performs worse than 
both proposed methods. For the case of different de¬ 
grees of noise (Fig. 3, top left) it can be seen that the 
total error increases with increasing noise. Similarly, in 
the case of different dimensions (Fig. 3, bottom left), 
the error increases with increasing dimensionality. This 
can be explained by the fact that the Frobenius norm in 
(23) sums over dP values. For various values of the graph 
density (Fig. 3, top right), the error for the H- and Z- 
matrix method is approximately constant (apart from 
the case of a rooted spanning tree at p = 0, according to 
Lemma 12.). 

5-4 Centralized methods for matrices in GL{d, M) 

In this set of experiments we compare the H-matrix 
method and the Z-matrix method. 

Using the reference-based method for the case of linear 
transformations is problematic because this method 
inverts the matrices Gij for {i,j) G £;mm-QSC ^ggg 
(23)). Therefore, for reasonably large noise, it is likely 
that there is some {i,j) € is ill- 

conditioned, resulting in the corresponding term in g' 
blowing up. In Fig. 5 this problem is illustrated, where 
the horizontal axis is shown in log-scale. The lines 
of the Z- and H-matrix methods almost coincide, so 
only the green line of the H-matrix method is visible. 
The reference-based method’s (black) line results in 
extremely large errors. Due to this reason, and since 
we have already shown that for the case of orthogonal 
transformations the reference-based method is inferior, 
in the following the reference based method is not used 
in the comparisons. 

For the complete graph case, in Figure 7 the improved 
performance the Gauss-Newton method, i.e., Algorithm 


n=30, d=3, p=0.5 



Fig. 5. Normalised error in (23) on the horizontal axis shown 
as log-scale for the Z-matrix method (blue), the H-ma- 
trix method (green) and the reference-based method (black) 
when considering transformations in GL{d, R). Note that the 
blue and green line coincide. 


a = 0.3, n = 10, p = 1 



Fig. 7. Performance of the Z-matrix method (blue), the 
i7-matrix method (green), and the Gauss-Newton method 
with the solution of the H-matrix method as initialization 
(black). 

3, (with the solution of the H-matrix method as initial¬ 
ization) is shown. The Gauss-Newton method run 5 iter¬ 
ations, but from inspection it could be deduced that the 
main convergence occurs already after two iterations. 

In Fig. 6, the comparisons of the H-matrix method and 
the Z-matrix method are shown. It can be seen that for 
small noise (Fig. 6, top left) both methods are compara¬ 
ble, whereas for a larger amount of noise the H-matrix 
method is able to obtain a smaller error. Similarly, for 
transformations with small dimensionality (Fig. 6, bot¬ 
tom left), both methods are comparable whereas for 
larger dimensions the gap between both approaches in¬ 
creases. On the contrary, (Fig. 6, top right) illustrates 
that with increasing graph density the line of the Z- 
matrix method approaches that of the H-matrix method 
(apart from the spanning tree case when p = 0, anal¬ 
ogous to the orthogonal transformation experiments). 
This indicates that the H-matrix method performs bet- 
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n=30, d=3, p=0.5 


a=0.3, n=30, d=3 



Fig. 3. Normalised error in (23) for the reference-based method (black), the Z-matrix method (blue) and the H-matrix method 
(green) when considering transformations in 0{d). In each sub-figure, a different parameter varies along the horizontal axis. 


n=30, d=3, p=0.5 


ct= 0.3, n=30, d=3 






Fig. 4. Gap function in (17) for the Z-matrix method (blue) and the H-matrix method (green) when considering transformations 
in 0(d). In each sub-figure, a different parameter varies along the horizontal axis. 


ter than the Z-matrix method if there is only little infor¬ 
mation available. A similar observation can be made for 
various n (Fig. 6, bottom right). For each subfigure, 100 
simulations for a certain configuration of tr, n, d, and p 
are shown. 


5.5 Methods for affine and Euclidean transformations 

In Figure 8 - for affine and Euclidean transformations 
- a comparison between four different methods can be 
found. The Gij transformations are affine respective Eu¬ 
clidean, but only the Algorithm 4 methods (red and 
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Fig. 6. Normalised error in (23) for the Z-matrix method (blue) and the H-matrix method (green) when considering transfor¬ 
mations in GL{d,'M.). In each sub-figure, a different parameter varies along the horizontal axis. 


<7 = 0.3, n = 10, p = 1 <7 = 0.3, n = 10, p = 1 




Fig. 8. Left figure: Performance of the methods for affine transformations. The Z-matrix method (blue), the iL-matrix method 
(green), the first two steps of Algorithm 4 (red), and Algorithm 4 (black). Right figure: Performance of the methods for 
Euclidean transformations. The Z-matrix method (blue), the //-matrix method (green), the first two steps of Algorithm 4 
where Algorithm 5 has been used instead of Algorithm 2 (red), and Algorithm 4 where Algorithm 5 has been used instead of 
Algorithm 2 (black). 

black) preserve this property. In the bottom right fig¬ 
ure the Gi transformations obtained in the Z-matrix 
method respective the //-matrix method have been pro¬ 
jected onto /?((/), i.e., the set of Euclidean transforma¬ 
tions. The orthogonal matrix part of the Gy transforma¬ 
tions were generated according to the description above. 

The elements in the transnational vectors were drawn 
from the uniform distribution over (—5,5) and addi¬ 
tional element-wise noise was added. 


5.6 Distributed methods 


Results of the distributed Z-matrix method are shown in 
Fig. 9 and results for the distributed H-matrix method 
are shown in Fig. 10. 
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Fig. 9. Normalised error in (23) on the vertical axis for the Z-matrix method (blue) and its distributed version (red) when 
considering transformations in 0{d). The horizontal axis shows the number of steps, where the step size has been chosen as 
e = 0.01 in each sub-figure. 


Conclusions 

This worked addressed transitive consistency of linear 
inverible transformations between Euclidean coordinate 
systems. Given a set of linear invertible transformations 
(or matrices) - that are not transitively consistent - 
the proposed methods synchronize the transformations. 
This means that they provide transformations that are 
both transitively consistent and close to the original 
non-synchronized transformations. First two different 
direct or centralized approaches were proposed. In 
the first approach - the Z-matrix approach - linear 
algebraic conditions were formulated that must hold 
for transitively consistent transformations. Then the 
sought transformations are obtained from the solution 
of a least squares problem. In the second approach - 
the i/-matrix approach - optimization problems were 
formulated directly, without taking a detour via linear 
algebraic constraints. The sought transformations are 
obtained from the solution of the optimization problems. 

A Gauss-Newton iterative method was also proposed 
where the solution from the iJ-matrix method was 
used as initialization. This method was later adapted 
to the case of affine and Euclidean transformations. It 


was shown in numerical simulations that for the case 
of affine and Euclidean transformations, this approach 
outperforms the i/-matrix approach and the Z-matrix 
approach. However, for orthogonal transformations no 
improvement is possible over the i7-matrix method. 

In a later part of the paper, for orthogonal matrices, two 
distributed algorithms were presented. These algorithms 
share similarities with linear consensus algorithms for 
distributed averaging. It was shown that these simple 
consensus-like protocols can be used to provide a solu¬ 
tion to our problem that is very close to the global opti¬ 
mum - even for noise large in magnitude. The proposed 
methods - both the direct/centralized and the itera¬ 
tive/distributed - have been verified to work in numer¬ 
ical experiments for a wide range of parameter settings. 
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